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Abstract 
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The formation of localized structures in the chlorine dioxide-idodine- 
Ch ' malonic acid (CDIMA) reaction-diffusion system is investigated numerically 

>^ . using a realistic model of this system. We analyze the one-dimensional pat- 

Q\ '. terns formed along the gradients imposed by boundary feeds, and study their 

linear stability to symmetry-breaking perturbations (Turing instability) in the 
plane transverse to these gradients. We establish that an often- invoked sim- 
^ ' pie local linear analysis which neglects longitudinal diffusion is inappropriate 

. for predicting the linear stability of these patterns. Using a fully nonuniform 

, analysis, we investigate the structure of the patterns formed along the gradi- 

I ents and their stability to transverse Turing pattern formation as a function 

00 ' of the values of two control parameters: the malonic acid feed concentration 

' and the size of the reactor in the dimension along the gradients. The results 

Q . from this investigation are compared with existing experiments. 
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I. INTRODUCTION 



Recent experimental developments have made possible the study of asymptotic spa- 
tiotemporal behavior in chemical systems in a controlled and reproducible manner, allowing 
predictions from theoretical and numerical studies of these systems to be compared quanti- 
tatively with experiments, in the same way that fluid systems have been studied. Indeed, 
the understanding of spatial pattern formation in nonequilibrium systems has greatly ben- 
efited from careful and controlled experiments on fluid systems [|I|. Unlike fluid systems 
which at high nonlinearity break down to a turbulent state characterized by a wide range 
of spatial scales, spatial patterns in chemical systems can be studied at high nonlinearity 
0], thus providing an opportunity to study rich and new phenomena that complement our 
knowledge from pattern formation in fluid systems. 

The symmetry-breaking instability of a system from a homogeneous state to a patterned 
state, predicted in 1952 by Alan Turing [^], was observed for the first time nearly 40 years 
later, in the chlorite-iodide-malonic acid (CIMA) reaction-diffusion system The Turing 

instability is characterized by an intrinsic wavelength resulting solely from reaction and 
diffusion processes. For this reason, it has particular relevance to pattern formation in 
biological systems 0. 

In contrast to hydrodynamic systems for which the governing equations and parameter 
values are well understood, how to model complex chemical systems is often not well known 
A realistic model of the simpler chlorine dioxide- iodine-malonic acid (CDIMA) reac- 
tion, which is similar to the CIMA reaction in terms of its stationary pattern-forming and 
dynamical behavior, has been proposed by Lengyel, Rabai and Epstein (henceforth referred 
to as LRE) P,^. Hence, the CDIMA reaction-diffusion system has the potential to become 
an archetype for the study of nonequilibrium pattern formation in chemical systems, 
in principle allowing numerical and theoretical investigations to be compared directly with 
experiments. 

In practice, however, this has not been fully realized for two reasons. First, numer- 
ical investigations of reaction-diffusion equations using realistic chemical parameters is a 
demanding computational task. In addition, the algebraic complexity of the realistic non- 
linear reaction terms renders these models unsuitable for analysis by standard analytical 
tools. As a result, the theoretical work on reaction-diffusion systems has been mostly based 
on abstract models. Second, despite the existence of the realistic CDIMA chemical model 
which has similar pattern-forming and dynamical properties to the related CIMA system, 
experimental work has continued to be based on the CIMA reaction, making direct compar- 
isons of numerical and analytical work with experiments difficult. Consequently, unlike in 
fluid systems, experimental and theoretical efforts in chemical systems have not been closely 
coupled. 

In this article, we use the realistic LRE model of the CDIMA reaction-diffusion system to 
investigate the Turing instability numerically |Tl| . Contrary to the case originally considered 



by Turing and subsequently by others, the experimental conditions under which Turing 
patterns form are not uniform, as required by the continuous feed of reservoir chemicals. We 
study the formation and stability of one-dimensional structures in the presence of boundary 
feed gradients. We first briefly review the Turing mechanism in Section 0. To facilitate 
comparisons with our numerical investigations, we describe the geometry and setup employed 
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by the relevant experiments in Section |IT|. The LRE chemical model is described in Section 
In Section |V A , we obtain the one- dimensional steady state chemical concentration 
profiles for a particular set of boundary conditions, and explore several different approaches 
to determine the linear stability of these profiles to transverse symmetry-breaking patterns. 
In Section |V lj| , the patterns along the gradients and their linear stability are further explored 
as a function of two control parameters. We summarize our results and consider prospects 
for further progress in Section |V1[ 



II. TURING INSTABILITY CONDITIONS 

In his original paper ||^, Turing suggested that the reaction and diffusion of chemicals 
could account for the instability of an originally homogeneous steady state to a stable steady 
pattern when triggered by random disturbances. These instability conditions, which are 
derived and discussed in detail in the text by Murray ^ are presented here for the purpose of 
introducing the notation for the rest of the article. Consider the general governing equations 
for the reaction-diffusion mechanism of two chemical species: 

a^ = fiu,v;fl) + V\ (1) 
— = g{u,v;fj.) + cV^v, (2) 

where / and g represent the (nonlinear) reaction kinetics, u and v are chemical concentra- 
tions, /2 is a set of reaction parameters that may include concentrations of other chemical 
species, c = D^/Du is the ratio of diffusion constants, and a > 1 is a constant separating the 
characteristic time scales for changes in the concentrations of the u and v species. Turing's 
idea was as follows if in the absence of diffusion {u{f, t), v{f, t)) tend to a linearly stable 
uniform steady state, then under certain conditions, the addition of diffusion leads to the 
development of spatially inhomogeneous patterns. Although these conditions were originally 
considered for a spatially uniform system, where the parameters fl are constant, the actual 
experimental realization of the Turing instability occurs in the presence of externally im- 
posed feed gradients, where fl = fl{z). In this section, we derive the Turing linear instability 
conditions for both uniform and nonuniform parameters, /I. 



A. Uniform background 

The uniform background case is realized experimentally in batch reactors where there 
are no externally imposed gradients from continuous feed of chemical reactants, and Turing 
patterns are necessarily transient. The parameters /7 are constants independent of position. 
The homogeneous steady state cl = {ug", Vs°) is obtained as the solution to: 

f{u,v;fl) = g{u,v;fl) = 0. (3) 

The linear stability of this state is obtained by substituting into the governing reaction- 
diffusion equations: 
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c(r,t) = Cs + Sc{r,t), 



(4) 
(5) 



where c = (m, v) is a vector of concentrations, k is the spatial wave number of the pertur- 
bation, Afc is the growth rate of the fc*^ mode, and {5u%, 5vl) is the corresponding constant 
eigenvector. The resulting linear eigenvalue problem: 



(an -k^) ja 012/0- 
022 - ck^ 



021 



6ul 



(6) 



yields, 
where 



2a 



[ac + \)k'^ - (an + 0-022) 



± 7^ V [(o^c + 1)P - (an + cra22)]' - 4/i(P), (7) 



/i(A;' 



a 



ck^ — (022 + cau)^ + (ana22 — «i2a2i' 



(8) 



The quantities an = /„, ai2 = fv, ^21 = Qu, and 022 = Qv are the elements of the Jacobian 
of the reaction terms with respect to the concentrations, evaluated at the uniform steady 
state. Afc has a rich behavior depending on the values of a, c, and aij. The conditions for the 
Turing instability are that this uniform steady state be linearly: (i) stable to homogeneous 
perturbation, and (ii) unstable to inhomogeneous perturbations. Hence, this is a symmetry- 
breaking mechanism, since it breaks the homogenous spatial symmetry of the uniform state. 
For the general reaction-diffusion system given in Eqs. (|l|) and (Q), these conditions are 
derived in Appendix ^ Below, we refer to the relevant results for the purpose of discussion. 
Stability of the uniform steady state to homogeneous k = perturbations requires the 
following inequalities to be satisfied: 



an + o-a22 < 0, 



(9) 



0110-22 - O12O21 > 0. 



(10) 



In order for the uniform steady state to be simultaneously unstable to inhomogeneous 7^ 
perturbations, we must have 



a22 + can > 0, 



(11) 



(a22 + can)^ — 4c(ana22 — 012021) > 0. 



(12) 

Comparing Eqs. (|^) and (|1T]) we conclude that an and 022 must have opposite sign. In 
the standard terminology, the activator species has a positive sign and the inhibitor has a 
negative sign in the Jacobian. Thus, taking an > and 022 < identifies u as the activator 
and V as the inhibitor. If cr = 1, then Eqs. (^ au' 
c > 1. In fact, given values of other parameters, c ^ 1 is required. 



TTD are simultaneously satisfied only for 
Since diffusion constants 
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of ions in aqueous solutions are all nearly the same ((9(10~^)cm^s~^), for the instability 
conditions to be satisfied, a must be greater than one. In Section the requirement 
(T > 1 will be put in the context of the fortuitous role of the starch color indicator in the 
pattern formation itself by providing the mechanism for slowing the activator reaction and 
diffusion with respect to those of the inhibitor. 

Equations @-(|l2D constitute the Turing conditions. It will be useful for future compar- 
ison of the local stability analysis with the full stability analysis of the nonuniform steady 
state to consider the Hopf bifurcation of the uniform system. For reaction parameters such 
that 

(an - 0-022)^ + 4crai2a2i < 0, (13) 

there will be a complex conjugate pair of eigenvalues for wave numbers in the range < 

fc^ < , where is given in Eq. ( |A12| ). With the above inequality satisfied, a Hopf 
bifurcation of the uniform system occurs when (an + (7022) > 0. Beyond the Hopf bifurcation 
point, there will be an unstable complex conjugate pair of eigenvalues for wave numbers in 
the range given by Eq. ( |A17| ). 



B. Nonuniform background 

In this case, the parameters /i, which depend on the concentrations of background chem- 
icals fed through the boundaries, are not constant but rather are functions of the variable 
z along the direction perpendicular to the feed boundaries. The steady state solution will 
now be a function of z, satisfying: 

f{usiz),Vsizy,fIiz)) + ^ = 0, (14) 

d'^Vs 

g{us{z),Vs{z);fI{z)) +c-^ = 0, (15) 

with Dirichlet boundary conditions at ^ = and z = L^. The stability of {us{z)^Vs{z)) is 
given by linearizing about this state: 

c{r,t) = Cs{z)+5c{r,t), (16) 
5c{r,t) = ^(5M,Jz),5t;,Jz))e*^^^-^^^eV*, (17) 

fc_L 

where fc^ is the wave vector perpendicular to the direction of the gradients. For simplicity 
of notation, we will drop the subscript "_L," taking k to be the transverse wave number. 
The resulting eigenvalue problem is: 



{aii{z) + ^ - k"^^ / a ai2{z)/a \ ( 5uk{z) \ f 5uk{z) 
a2i{z) a22{z) +c£-ck' [ Sv^iz) J " [ Sv^iz) 



(18) 



with {5uk{z),Svk{z)) satisfying the same Dirichlet boundary conditions as the steady state. 
As Pearson et al. have noted, this is an infinite-dimensional eigenvalue problem for 



5 



each k which is formally similar to the Schrodinger equation. However, the Jacobian of 
the reaction terms is not symmetric, rendering the linear operator non-Hermitian. Hence, 
it must be solved numerically, by discretizing the z spatial direction into A^^ niesh points 
and solving the resulting 2Nz x 2Nz matrix eigenvalue problem for each k. This method of 
solution is described in Section |V A3. 



C. Locally uniform background 

In the presence of ramps in control parameters, a naive assumption is that a structure 
will form in the region of space where the local value of the control parameter allows it to be 



stable in the corresponding uniform problem |15|. This "locally uniform" approach amounts 



to treating each location along the gradients in the z-direction to be an independent and 
uniform quasi-two-dimensional system in the x — y plane. The corresponding locally uniform 
steady state which depends parametrically on z is given by the solution to: 

f{u,v;fl{z)) = g{u,v;fl{z)) = 0. (19) 

The Turing instability conditions can then be examined at each point in z to determine 
whether or not a linear analysis predicts the formation of transverse Turing patterns in any 
interval along z. 

Since the resulting eigenvalue problem for the stability of the locally homogeneous steady 
state to a symmetry-breaking instability requires only a 2 x 2 analysis at each z, it is 
computationally simple. The validity of this local analysis is assessed in Sec. |V A2|, by 



comparing the result with that from the fully nonuniform analysis (Eq. |T^) of the steady 
state along the gradients. 

III. EXPERIMENTAL GEOMETRY 

The first experimental realization of the steady-state Turing patterns predicted in 1952 
was made in 1990 by Castets et al. [Q, and was subsequently confirmed by others 
This was made possible by the development of open spatial reactors which allowed ex- 
perimentalists to maintain a reaction far from equilibrium through a continuous supply of 
reactants, while avoiding convective transport. These sustained patterns have been obtained 
in only one controlled experimental system to date, the chlorite-iodide-malonic acid (CIMA) 
chemical react ion- diffusion system. The principles of operation of these reactors have been 
discussed elsewhere in detail @,|,|^,|l^ . 



In this section, we introduce the thin-strip reactor which is investigated numerically 
in this work, since a geometrical description of the experimental setup is useful in the 
understanding of our results. A detailed description of the chemical model is presented later 
in Section and is not necessary for the discussion presented here. The thin strip reactor 
is comprised of a thin rectangular gel strip, such that L ^ w > /z, as in Figure |1]. Typically, 
h <1 mm, L ~ 25 mm, and tf ~ 3 mm. The gels are water-based, acting as essentially water 
in a loose polymer grid. The gel core of the reactor is in contact with two continuously stirred 
reservoirs of chemicals. Components of the reaction are distributed in the two reservoirs in 
such a way that neither is separately reactive. In the CIMA experiments, these reservoir 
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species are malonic acid (CH2(COOH)2 or MA), iodide and chlorite (C102~). The LRE 
model of the CDIMA system takes as input malonic acid, iodine (I2) and chlorine dioxide 
(CIO2). As the reservoir species diffuse and react through the gel, the two dynamical species, 
iodide and chlorite, which take part in the pattern-forming instability, are produced. The gel 
is preloaded with a soluble starch that acts as an indicator by changing color from yellow to 
purple with changes in triiodide concentration. The large starch molecules are immobilized 
in the gel matrix, and for this reason actually play a role in the pattern formation itself. 

Observations are made in the direction perpendicular to the x—z plane (along Oy), which 
allows viewing of the mulitfront patterns that develop along the boundary feed gradients (in 
the z-direction), as well as patterns that form parallel to the boundaries (in the a;-direction), 
breaking the boundary feed symmetry (single or multiple layers of spots). The symmetry- 
breaking patterns form in a thickness A along z. If the gel is thin enough (/i ~ A of the 
Turing patterns), these patterns are one- or two-dimensional, depending on whether A is 
of the order of one or more wavelengths A of the structure. With h ^ X, for example 
/i ~ L as is the case with disc reactors, the patterns are quasi-two-dimensional (referred to 
as "monolayers") for A ~ A, or three-dimensional for A > A (referred to as "bilayers" for 
A ~ 2A). 

A modified thin-strip reactor, where the feed surfaces {L x h) are no longer parallel but 
make an angle, has been developed and used by Dulos et al. |^^, where h = 0.2 mm, L = 25 
mm, and the w ranges from 1.75 to 3.5 mm. The variation in w causes a gradual change in 
the reservoir concentration ramps across the gel, the effect of which can be studied on the 
patterns that form along and transverse to the gradients. 



IV. CHEMICAL MODEL OF THE CDIMA SYSTEM 



Lengyel, Epstein and Rabai have proposed a model for the temporal oscillations in the 
chlorite-iodide-malonic acid reaction, C102^-I~-MA , which is based on the simpler chlorine 
dioxide-iodine-malonic acid reaction, CIO2-I2-MA , referred to as CDIMA P,^. They have 
shown experimentally that the CDIMA system also exhibits the Turing instability in both 
closed and open systems |T^J21|]. By monitoring the CIMA reaction in a closed system 
spectrophotometrically, it was determined that after an initial fast consumption of I~ and 
C102~ during a pre-oscillatory period to produce I2 and CIO2, the reaction of CIO2, I2 and 
MA accounts for the oscillations. The LRE CDIMA model consists of three reactions for 
MA, I2, CIO2, I~, C102~ and |T^, with empirically determined rate laws: 



MA + I2 ^ IMA + r 
d[h] _ h,[MA][h] 
dt 

CIO2 + r 

_rf[C102] 

dt 



ri 



kib + [h] 

^ CIO2" + 1/2 I2 
= A;2[C102][r] = r2 



(20) 



(21) 



CIO2" + 4r + 4H+ ^ cr + 2I2 + 2H2O 
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rft ' ' h+[l- 



Lengyel et al. ||18[ have modeled the effect of unreactive starch-complex formation on 
the CDIMA system, where the complexing agent is (S + I2). Although formation of the 
starch-triiodide complex (Sla") is a complicated process, it can nevertheless be described as 
a single reaction: 

S + I- + I.^Sl3-. A-^^^i±. (23) 

where K is the equilibrium constant, and the reaction rate is given by: 

r4 = A;+[S][l2][r]-A:_[Sl3-]. (24) 

Using the above, the full reaction-diffusion model for the CDIMA system, with the addition 
of the reaction with starch, is given by ^ 

d[MA] 



dt ^ 2 

dt 



dt 

d[C\02~ 

dt 
d[YL+ 



ri + Z}MAV'[MA] (25) 

n + + 2r3 - r4 + D,, [I2] (26) 

r2 + /^cio.V2[C102] (27) 

ri - r2 - 4r3 - r4 + A- [!"] (28) 

r-2-r3 + Z^cio.-V2[C102l (29) 

n (30) 

ri-4r3 + Z^H+V2[H+]. (31) 

The rate and diffusion constants used in the numerical calculations here are taken from Refs. 
|T3,|T|,g§ and are given in Table |. 

Lengyel et al. Q have shown that these reactions successfully simulate the temporal 
behavior of [MA], [I2], [CIO2], and [C102~] in a batch experimental system. Their 
numerical results show that while the intermediates and [C102~] vary by several orders 
of magnitude during an oscillation, [MA], [I2] and [CIO2] vary more slowly. In addition, the 
contribution of [H+] to the reaction terms is relatively small, and this species can be ne- 
glected. This suggests a reduction of the full model to a three- variable system ([I~], [C102~] 
and [Sis"]) by treating the concentrations of the [MA], [I2] and [CIO2] reactants as con- 
stants, making the model mathematically and numerically more tractable. This procedure 
illustrates the adiabatic elimination of fast modes in dynamical systems which reduces the 
full dynamics to only a few degrees of freedom. 

It is further assumed that there is a large excess of starch uniformly distributed so that 
its concentration is always very close to its initial value [S]^, and that the complex formation 
and dissociation is fast. Then, 
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[Sl3-]^i^[S][y-[I-]^i^'[I-], (32) 

where K' = iir[S]Q[l2]o. Adding Eqs. ( p8|) and (^), and using Eq. (|3^) , a two-component 
reaction-diffusion system is obtained: 

0-— = fci -k2U- — — - + DuV u (33) 
ot h + 

= k2U- — — - + D^V V, (34) 
ot h + 

where fc/ = A;i4MA]4l2]o/(fcif, + [l2]o), k2' = A;2[C102]o, W = k^blhU and a = 1 + K' > 1. 
The subscript "o" refers to the concentrations of species which are taken to be constant, 
and u and v represent the concentrations of I~ and C102~ species. The role of the immobile 
starch color indicator in providing the relative slowdown of the reaction and diffusion of the 
activator with respect to that of the inhibitor enters through the parameter a > 1. 



V. ANALYSIS 

In this section, we use the LRE chemical model to investigate several aspects of the 
experimental CIMA system. The focus is to demonstrate the potential for quantitative 
analysis of experimental results using the realistic CDIMA chemical model. In Section |V A|, 



we investigate numerically the formation of one-dimensional multifront localized structures 
along the gradients of imposed boundary feeds. We study the linear stability of these 
structures to transverse symmetry-breaking perturbations using the two-variable reduction 
of the LRE model. We compare our results from a local analysis to that from a fully 
nonuniform analysis. We review a proposed modification to the local analysis and show that 
it does not successfully account for the presence of gradients. In Section |VB| , we further 
explore the structure and linear stability of the one- dimensional patterns along the boundary 
feeds as a function of two control parameters: the malonic acid reservoir concentration and 
gel width. We map out the linear instability intervals in each case. We discuss the qualitative 
agreement of our results with relevant experiments. 

A. Linear Analysis of One-Dimensional Patterns Along Gradients 

1. Stationary solution along the z-direction 

The full 7-component LRE model equations given in Eqs. (|25|)-(pID were evolved forward 
in time to obtain the steady state solution in one dimension along the gradients. The 
boundary conditions, [MA]^^ = 1 x 10"^ M at the left boundary, and [12]^ = 1 x 10"^ M 
and [0102]^ = 6 X 10~^ M at the right boundary, were chosen so as to be consistent with a 



previous numerical investigation of the LRE model in one dimension by Lengyel et al. [jT9 
Since the boundary conditions giving a transverse instability were a priori unknown, we used 
these values as our starting point. The spatial z-direction was discretized on an irregularly 
spaced mesh to allow a greater number of mesh points in the regions where there was more 
structure in the solution, without excessively increasing the overall number of mesh points 
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in the problem. A five-point finite- difference approximation to tlie diffusion operator was 
used on tlie variable mesh. The numerical scheme employed for the time evolution was 
Crank-Nicholson implicit time-stepping for the linear terms, and Adams-Bashford explicit 
time-stepping for the nonlinear terms. A banded solver ||2^ was used at each time step to 
solve for the solution at the next time step. The initial concentrations were uniform in the 
z-direction (and equal to 5 x 10^^^ M). The time evolution was continued until there was 
no appreciable change in the solution. 

The results are displayed in Figure ^ The steady state solution for the starch-triiodide 
complex (Sis") plotted in Figure |^(f) represents the experimentally observable profile. As 
expected, it tracks the iodide (I~) steady state solution in Figure 0(d). It corresponds to 
a (non-symmetry breaking) pattern of stripes parallel to the feed boundaries (in the x — z 



plane, with x being the uniform direction), such as those observed by Perraud et al. [jT5 
although the boundary species are different in their experiments on the CIMA reaction from 
those considered here. Upon increasing the [C102~] reservoir concentration, they observe a 
break-up of the stripes to rows of spots parallel to the feed surfaces. This is a symmetry- 
breaking instability, since the boundary feed symmetry of the system is broken. In the 
following, we will investigate the linear stability of our numerical steady state along the 
gradients to such a transverse pattern-forming instability. 



2. Locally uniform stability analysis 

To examine the stability of the stationary patterns that form along the gradients of 
boundary feeds (z-direction), the simplest approach is to treat each location z as being 
independent and locally uniform in the transverse plane (see Section p^I C| ), thereby neglecting 
diffusion along the 2;-direction (longitudinal diffusion). The locally homogeneous stationary 
state in the dynamical species, I~ and C102~, at each z can be constructed either from the 
linear (diffusion only) profiles of the reservoir species, MA, I2, and CIO2, or more correctly. 



from their reaction-|-diffusion profiles, obtained by evolving the full model, Eq. (25)-(pTD. 
In either case, using the two-variable activator-inhibitor reduction of the LRE model, Eqs. 
(p3|)-(|3^), the resulting eigenvalue problem for the stability of the locally homogeneous 
steady state to a symmetry-breaking instability requires a simple 2x2 analysis. Hence, 
it is desirable to use such a local approach if it can be shown that it accurately describes 
the physical problem. In that transverse instability would occur in a region of 

width A along the gradients which is linearly unstable to a A; 7^ instability. Indeed, even 
though Turing patterns are obtained under experimental conditions that by necessity lead to 
nonhomogeneous parameter ramps, a local linear analysis is most commonly used to predict 
the formation of a transverse symmetry-breaking instability. In this section, we examine in 
the context of the two- variable LRE model the locally uniform approach to determining the 
stability of the stationary patterns that form along the gradients of boundary feeds. 

The locally uniform steady state in the variables I~ and C102~ at each point in z along 
the gradients of the background chemicals is shown in Figure 0. This solution is obtained 
according to Eqs. (p3D -(^) using the numerical reaction+diffusion profiles of Figure ^ for 
the MA, I2 and CIO2 species, but neglecting the diffusion terms. The dependence on z in 
this plot is parametric. 

The stability of the local steady state at each z is obtained from Eq. (^. This analysis 
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predicts the existence of a finite instability region. Tlie curves in Figure ^ represent tlie 
Turing instability condition boundaries, Eqs. (|9|)- (p!2D , for each locally uniform steady state. 
We have also plotted the boundary which, when less than zero, gives the interval in z where 
the locally uniform steady state has a complex conjugate pair of eigenvalues for a range in 
k given by Eq. ( |A12| ). The vertical axis has no labels since we have plotted quantities which 
have different dimensions and only their signs are of interest. The shaded region denotes 
the interval over which the locally uniform steady state in the x — y plane is linearly stable 
to homogeneous perturbations and unstable to inhomogeneous perturbations. The width 
A of this region is approximately 0.15 mm, which is within the 0.13 — 0.33 mm range of 
experimentally observed Turing wavelengths p3| . 

Figure ^ shows the gain curves for the locally uniform steady states at selected points 
along the z-axis, consistent with the above linear stability boundaries. Note that Figs. 
^(d)-(i) illustrate the role of the complexing agent (S + 12) in suppressing the oscillatory 
instability, since the concentration of I2 (and therefore the complexing strength) sharply 
increases as the right boundary is approached. 



3. Fully nonuniform stability analysis 

To assess the validity of the locally uniform stability analysis presented above, we have 
carried out a fully nonuniform analysis, as described in Section |1IB| . The eigenvalue problem 
given in Eq. (|I8|) was discretized on the same variable mesh as that on which the nonuniform 
steady state was obtained, and was solved for all eigenvalues and eigenvectors at each value 
of transverse wave number k using EISPACK [^|. Since we have a general real matrix. 



with no special features such as symmetry, the most general routines were used. Figure 
^ shows the real part of successive eigenvalues with largest real parts. This result shows 
the steady state to be stable to all transverse perturbations. The eigenvector with slowest 
decaying (real) growth rate at k = 81.6 cm~^ is plotted in Figure 0. It is locahzed roughly 
in the region along z where the locally uniform stability analysis predicts the corresponding 
uniform steady state to be unstable to a transverse instability. 

Since we are generally interested in the most unstable mode, in this case, we checked the 
numerical validity of the gain curve for the slowest decaying eigenvector (top-most continuous 
curve in Figure ^ against two different numerical methods. First, the linear system for 2Nz 
variables (eigenvector), where A^^ is the number of mesh points, was solved as a nonlinear 
root finding problem in {2Nz + 1) variables, including the eigenvalue. Second, starting with 
the eigenvalue and eigenvector based on the previous two methods, inverse iteration was 
used to verify the results. Both checks agree with the results from EISPACK. 

The convergence of the slowest decaying gain curve as a function of the variable mesh 
was also investigated. Starting with a particular distribution of mesh points, the mesh size 
was successively halved, the corresponding steady state obtained, and the eigenvalues and 
eigenvectors found. Basically, the number and distribution of mesh points must be sufficient 
to well resolve both the structure of the steady state and its most unstable eigenvectors 
for numerical convergence. All numerical calculations were performed on an IBM RS6000 
workstation, with the exception of eigenvalue/eigenvector determination using EISPACK with 
greater than approximately 500 mesh points, which was done on a CRAY C90. 
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This result contradicts that from the locally uniform analysis which predicts a linear 
instability for these reaction parameters and boundary conditions. Below, we review a 
proposed modification to the locally uniform analysis, and determine whether it is sufficient 
to bring the local analysis closer to the fully nonuniform one. 

In the above analysis and in those discussed in Section |V A 2| , we have used the two- 
variable reduction of the LRE model. We have directly verified the two-variable reduction 
of the full 7- variable (including H"*") model by performing a 7-variable linear stability analysis 
of the nonuniform stationary state using inverse iteration. By comparing the 7-variable and 
2- variable linear stability results, we have implicitly verified the assumption that the reservoir 
species, MA, I2 and CIO2, do not play a role in determining the pattern-forming instability 



of the stationary states [11 



4. Modified local stability analysis 

The locally uniform analysis neglects diffusion along z which couples quasi-two- 
dimensional uniform slices. Lengyel, Kadar and Epstein (LKE) have proposed a mod- 
ification to account for this diffusion, assuming that diffusion along z is relevant only on a 
length scale of the order of the Turing wavelength. The basic idea behind the LKE modifi- 
cation is simple. In the presence of gradients in the 2;-direction, the Turing unstable mode 
is "split" between its "longitudinal" (along z) and "transverse" dependence: 

= + ^i.) (35) 

where kc is the critical wave number in the (narrow) Turing unstable region along z, depend- 
ing only on the local values of reaction and diffusion parameters. A transverse instability 
can occur provided the width of this region is not smaller than a Turing wavelength. 

This modified local analysis is used to better predict the region along the gradients where 
a transverse instability occurs, and to obtain more accurately parameter values for inves- 
tigating (transient) Turing patterns in batch reactors. The mechanics of the modification 
consists of adding an approximate term to the governing equations for the diffusion of the 
steady state along z, which does not alter the composition of locally uniform steady state 
but does affect its stability. This approximation to the diffusion operator is given by: 

d^u u{z - ^) - 2u{z) + u{z + ^) 8{u-u{z)) 

(f)^ ' ^ ^ 

u is the average value of the locally uniform steady state on the two sides of the region 
of width A, which is characteristic of the longitudinal variation of the steady state. The 
validity of this estimate relies on the smallness of this width. The reaction terms / and g 
are modified: 

f{u, V- fl{z)) = f{u, V- fl{z)) + 8D^{{u - u{z))/A^ (37) 
g'{u, V- pi{z)) = g{u, v; pi{z)) + 8D,{{v - v{z))//\\ (38) 

and the Jacobian of the reaction terms, Ojj, in the linear stability analysis is modified 
accordingly. The Turing conditions can be rewritten as: 
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where 



K'{z)>Hi{z)>H2{z), 



(39) 



Hi = -aii/a22 - 1, (40) 
On 

2wc(aiia22 — 012^21) — can 



^ , 1. (41) 



The range of the Turing instabihty is given by the crosspoints of K' and Hi and of Hi and 
H2. Since aij depend on A which is a priori unknown, Hi and H2 are evaluated iteratively 
from an initial estimate for A until convergence is achieved. 

Figure |^ shows that the effect of this modification is to extend the z-range of the trans- 
verse instability from approximately one to two Turing wavelengths. The boundaries given 
by the functions K'{z), Hi{z) and H2{z) are combinations of the boundaries given in Figure 
^ which resulted directly from the linear instability conditions. Therefore, although the 
representation of the boundaries in Figure ^ differs from that in Figure the instability 
region is the same in the unmodified case. 

The proposed modification increases the range of the Turing instability by suppressing 
the homogeneous oscillatory instability (moving the left boundary to the left), while not 
affecting the right boundary corresponding to the inhomogeneous instability. Identifying the 
left boundary of the 2;-range in which transverse patterns would form with the homogeneous 
instability is unphysical, since there is no mixing of modes at the linear level. Instead, it 
seems more appropriate to identify the left boundary of the Turing region with the 2;-location 
where an inhomogeneous instability ceases to exist (can + a22 = 0). 

An alternate modification to the local analysis is to carry out the linear stability anal- 
ysis about the nonuniform steady state along z, such that the diffusion operator in the 
governing equations acts only on the steady state and not on the instability eigenvector in 
the z-direction. This "local" analysis incorporates the effect of diffusion along z through 
the nonuniform steady state, but the instability eigenvectors are "local" and depend only 
parametrically on z. The resulting eigenvalue problem becomes: 

/ {aii{z) - e) /a ai2{z)/a \ ( 6uk{z) \ . Suk{z) \ 

V «2i(^) a22{z) -ck^ J \6vk{z) J ''^ ' \6vk{z) J ' ^ ' 

where aij{z) are evaluated at the nonuniform steady state, and {5uk{z) , 5vk{z)) and \k{z) 
depend parametrically on z. The stability boundaries are very irregular and not shown 
in this case. Except at the sharp edges of the nonuniform steady state and over a region 
roughly equal to the width of the sharp edge (much smaller than a Turing wavelength) , this 
analysis predicts no 7^ instability. 



5. Discussion 

We conclude that to accurately predict the formation and location of the Turing insta- 
bility region, at least the one-dimensional steady state along the gradients must be solved 
for numerically using the full model including longitudinal diffusion. A "local" stability 
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analysis about this steady state does reproduce the result from the fully nonuniform sta- 
bility analysis. However, a local analysis neglecting longitudinal diffusion of the stationary 
state does not correctly describe the linear stability of this state. We have presented here 
a first direct demonstration of this point by carrying out a fully nonuniform as well as a 
local linear stability analysis. As has been suggested two-dimensional (nonlinear) time 
evolution of the model is the definitive method for predicting a transverse instability. We 
have accomplished this for the LRE model, and the results will be published elsewhere ^9 



The locally uniform steady state profiles for the two dynamical variables I~ and C102~ 
(Figure ^) do not include diffusion along the z direction and are qualitatively different from 
the numerical solution including diffusion (Figs. 0(d), (e)). Hence, it is not surprising that 
the local stability analysis about this steady state does not agree with the fully nonuniform 
one. In particular, at the left boundary, the locally uniform C102~ profile is several orders 
of magnitude greater than the corresponding numerical solution including diffusion. This 
large discrepancy is accounted for by the diffusion of this species in a region extending over 
approximately the left half of the gel, as can be seen from the almost linear (diffusive) profile 
for C102~ over this region (Figure 0(e)). The LKE modification to the locally uniform 
analysis, which corrects for this diffusion of the steady state along z, assumes that it is 
relevant only on the length scale of the order of the Turing wavelength. For the parameter 
values investigated here, this assumption is not valid, and the modified stability analysis 
does not correctly predict the stability of the structures along the gradients. 

The dependence of the locally uniform steady state profiles for the intermediate species on 
the background concentration gradients is through the reaction terms. The results presented 
here are for background profiles obtained from reaction and diffusion of all six species. We 
have checked that these background stationary profiles do not vary considerably from those 
obtained by setting the intermediate species identically equal to zero. In this way, we rule 
out the possibility that the diffusion of the intermediate species feeds back into the profiles of 
the background variables, thereby accounting for the large discrepancy between the locally 
uniform steady state profiles of the intermediate species and those including diffusion. In 
particular, the large left boundary value of the locally uniform C102~ species results from 
the strong suppression of I2 relative to CIO2 at this boundary, which becomes even more 
pronounced with backgrounds obtained from the intermediate species set identically equal 
to zero. 

It is desirable to obtain semi-analytical solutions to the stationary structures along the 
gradients, which could then be used in (semi-analytical) linear stability analyses of these 
states. This has been done, for example, for the Brusselator model in the presence of 
slow spatial gradients using a WKB-like approach [^,^|. The localized structures along 
the gradients are obtained as marginally stable perturbations to the locally uniform steady 
state. However, in the case presented here, it is not possible to carry out a similar analysis. 
First, there is a large discrepancy at the boundaries between the locally uniform solution 
and the desired solution including diffusion satisfying Dirichlet boundary conditions. Even 
if the boundaries are ignored and an approximate solution in the interior of the gel is sought, 
our numerical results show that the steady state including diffusion is not a weakly nonlinear 
perturbation to the locally uniform steady state. Therefore, seeking a correction given by 
marginally unstable modes is not justified. For small gel width (see Sec. [V B 2|) , where the 
background concentration profiles are almost linear, such a WKB-like description relying on 
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the slowness of the parameter ramps can perhaps be sought. 

Numerical calculations based on the two-variable LRE model with uniform background 
have shown the transition to a symmetry-breaking instability to be strongly subcritical 



27| , |28[| . Although it is not clear how the range of parameters investigated in these works 
compares with their local values in the actual ramped experimental system or in our numer- 
ical example, these results imply that a linear stability analysis of the locally homogeneous 
steady states would not predict the existence of a finite amplitude instability in the sub- 
critical regime. The nature of the transition of the fully nonuniform stationary structures 
along the ramps to a transverse symmetry-breaking instability has not been determined 
yet. Should this transition be supercritical, or even weakly subcritical, the fully nonuniform 
linear stability analysis would well predict the formation of quasi-one- dimensional symmetry- 
breaking spots in the thin-strip experiments. This is currently under investigation [E3 . 



B. Parameter Dependence of the Turing Instability in the CDIMA System 

As discussed in the previous section, by adopting a locally-uniform approach, existing 
numerical studies of the stability of the stationary patterns along feed gradients have not 
fully taken into account the reaction+diffusion feed gradients. Thus, the parameter range 
for the occurrence of a transverse symmetry-breaking instability in a gradient system within 
the context of the LRE model is essentially unknown. Hence, we have undertaken a sys- 
tematic search, using the concentration of one of the reservoir species and the gel width as 
the control parameters. Variation of either of the control parameters changes the nonlinear 
reaction+diffusion profile of the background species, however they are not equivalent opera- 
tions. In the following sections, we present our numerical results and make connection with 
relevant experimental work. 



1. Variable malonic acid boundary condition 

The parameter search for the Turing instability in the CDIMA react ion- diffusion system 
as a function of the malonic acid concentration at the left boundary was performed for 
[MA]l in the range 0.004 M to 0.035 M. The concentrations [h]R and [ClOa]/? at the right 
boundary were held fixed at 0.008 M and 0.006 M, respectively. These values were chosen so 
as to lie within the range of the initial concentrations of these species used in experiments 
on this system in batch reactors [O , and therefore should also be experimentally accessible 



in open reactors. 

First, we numerically obtained steady state solutions of the full 7- variable governing 
equations as a function of [MA] 2,, as described in Section |V A 1| . The analysis described 
in Section |V A3| of the linear stability along the gradients to transverse symmetry-breaking 
perturbations was repeated for each stationary state. This was performed using the reduction 
of the full LRE model to the two dynamical variables I~ and C102~. The eigenvalue and 
eigenvector corresponding to the fastest growing (or slowest decaying) mode at each value 
of the transverse wave number k were obtained numerically using inverse iteration, and 
confirmed for select values of k using EISPACK ^ . 
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In Figure we have plotted the value of the control parameter [MA] 2, versus wave 
number, with the solid boundaries denoting marginally stable wave numbers. The shad- 
ing indicates the Turing-unstable regions. We note that the unstable regions are disjoint, 
corresponding to the following scenario: as the control parameter is continuously varied, 
the stable stationary state along the gradients first becomes unstable to transverse Turing 
patterns at a critical value of the control parameter, and initially remains unstable as the 
control parameter continues to increase. It becomes stable again once the control parameter 
exceeds a second and higher critical value. This is qualitatively consistent with the exper- 
imental observations of Perraud et al. []T3| on the CIMA react ion- diffusion system. Their 
results show that as the concentration of [CIO2 ]_r at the right boundary is increased, the 
number of alternating dark and bright bands parallel to the feed boundaries increases, and 
several layers break up into rows of spots. As [C102~]_r continues to be increased, the spot 
patterns develop along more bright bands, until they eventually disappear and the parallel 
stripes are recovered. 

In Figure |T0|, we have plotted the stationary solutions for the experimentally observed 
SIs^ species from our numerical calculation, at various values of [MA]/,. We observe that 
the number of peaks in the steady state solution increases with increasing [MA]^, while the 
characteristic "wavelength" of the patterns along the gradients remains relatively constant 
(excluding the leftmost peak). The transition from one unstable region to another corre- 
sponds roughly to the appearance of an additional peak in the stationary solution along 
the gradient: instability region I corresponds to a steady state with three peaks, region II 
corresponds to four peaks, and region III corresponds to five peaks. As an example, for the 
instability region II, we show in Figure |11| density plots of the steady state and the fastest 
growing (slowest decaying) wave vector. This eigenvector is localized at and roughly tracks 
one of the minima in the steady state solution. This is also the case for regions I and III. 
Specifically, the unstable vector appears to be approximately localized at the next-to-last 
minimum of the Sla" solution. The appearance of the symmetry-breaking instability at a 
minimum of the starch-triiodide is consistent with the above experimental observations. 

To better quantify this trend, in Figure 12, we plot the value of all six chemical species 
at successive minima of the Sla" stationary solution as a function of [MA]^. In this figure, 
the circles, triangles, squares and diamonds correspond to the second, third, fourth and 
fifth minima, respectively. The points corresponding to critical [MA] 2, values are filled (two 
points for each unstable region). We note that the relevant chemical species for tracking 
the instability are the two dynamical species, I~ and C102~, plotted in Figure |T^(d) and 
(e). (The Sis" concentration depends on the product of the I~ and I2 concentrations.) The 
concentration of I~ remains within the range of approximately (1.0 — 2.2) x lO^^M, and the 
concentration of C102~ remains approximately constant at 1 x 10~® in each of the instability 
regions. (Note that the concentration of Sis" does not stay within a more-or-less constant 
range for each of the instability regions due to the drift in the I2 concentration.) A simple 
interpretation of these numerical results is that an instability occurs when the concentrations 
of the dynamical species I~ and C102~ lie within a certain range. As the stationary solution 
changes with increasing [MA]/:, the instability vanishes when the value of the concentration 
of the steady state falls outside of this range, and reappears when the concentration at the 
next minimum is within this range again. 

The experiments of Perraud et al. [ITH] show that as the control parameter is increased, 
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the stationary pattern of stripes along the gradients first becomes unstable in a single stripe 
region and then in multiple stripes as the control parameter continues to increase. In our 
numerical investigations, the most unstable vector remains singly peaked in all cases. To 
further investigate this point, we have examined the spatial structure (along z) of the linear 
instability eigenvectors as a function of k. Figure |T3| shows a density plot of the most unstable 
(or least stable) eigenvector, corresponding to the I~ species, as a function of transverse wave 
number k for [MA]^^ = 0.023 M. The horizontal axis represents the spatial coordinate along 
the gradients, and the vertical axis is the transverse wave number, ranging from k = cm~^ 
to k = 900 cm-^ 

In Figure |T^, we show the first fourteen eigenvalues with largest real parts as a function of 
k. We note the eigenvalue crossings that define distinct "modes" cutting across the spectrum 
of eigenvalues. As an interesting aside, we have investigated the eigenvectors corresponding 
to one such "mode", denoted by filled triangles in this figure. The spatial profile of these 
eigenvectors changes continuously with increasing k. Hence, their interpretation as a single 
physical "mode" is not obvious. There does not appear to be a mode crossing between the 
first and second branches of Figure although the structure of the first eigenvector changes 
considerably from being multipeaked (and stable) to being singly peaked (and unstable, 
but subsequently stable again) as k increases. The experiments of Perraud et al. could 
correspond to a case where, as the control parameter is varied, a mode with multiple peaks 
becomes unstable. This does not occur in our numerical investigation, where the multiply 
peaked modes remain stable. The appearance of multiple instability layers could also be a 
nonlinear effect, resulting from linear growth and nonlinear saturation of the singly peaked 
unstable eigenvector, as well as growth of its side lobes. However, we have shown that 
in the two-dimensional evolution of the governing equations, the spatial profile of the most 
linearly unstable mode is in fact preserved in the LRE model. 



2. Variable reservoir length along boundary feeds 



The experiments of Dulos et al. [|14| have aimed at elucidating the transition from quasi- 
two-dimensional to three-dimensional Turing patterns by combining observations in bevelled 
thin-strip and disc reactors. Motivated by these experiments, we have undertaken a simi- 
lar numerical investigation which does not directly address the same question, but rather 
continues to focus on the localized patterns along the ramps. In particular, we consider the 
experimental results from the variable-width thin-strip reactor. In these experiments, the 
transition between the domains with one and two rows of spots, and the possible infiuence 
of the feed gradients on the phase relations between the spots in the two rows have been 
studied. The vanishing amplitude of the spot modulations before this splitting occurs is not 
well understood. 

Our numerical results address the latter question. We have reproduced the observed 
qualitative trend of the symmetry-breaking instability occurring and subsequently disap- 
pearing in a single layer as the gel width is varied. The boundary conditions are held fixed 
at [MA\l = 0.023 M, [hU = 0.008 M and [ClOajij = 0.006 M, while the gel width is varied 
from 0.14 — 0.39 cm. The stationary localized patterns along the gradients as a function of 
the scaled gel width are shown in Figure |T5|. Figure |1^ shows the linear stability of each 
solution to a symmetry-breaking instability. In Figure |T^, we have plotted the value of the 
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control parameter w versus wave number k, with the sohd boundary denoting the marginally 
stable wave numbers. For the parameter values and boundary conditions numerically inves- 
tigated here, we do not observe the transition from an unstable monolayer to an unstable 
bilayer as the gel width is increased. We have followed the most linearly unstable k ^ 
mode as the gel width is varied, and it remains singly peaked. 



These numerical results are analogous to those presented in Section V_B1, where the 



parameter dependence of the Turing instability in a ramped system was investigated as a 
function of the malonic acid boundary condition. Here also, we track the concentrations of 
the chemical species at the location of the minimum of the starch-triiodide complex where 
the fastest growing/slowest decaying instability eigenvector is localized. The results are 
given in Figure |18|: as before, triangles, circles, squares, and diamonds represent the second, 
third, fourth and fifth minimum, respectively. The solid squares correspond to the critical 
lengths, below and above which the instability vanishes. We note the following about the 
concentrations of the chemical species in the unstable interval (at the third minimum of 
the stationary SIs^ solution): (1) the concentration of C102~ is approximately at 1 x 10~^ 
M, in agreement with the variable [MA] 2, investigation, (2) the concentration of I~ is in 
the approximate range of 1.8 — 2.2 x 10~^ M, again in agreement with the variable [MA]^ 
investigation, and (3) the concentrations of the background species, MA, I2 and CIO2 are 
approximately equal to those in the variable [MA] 2, case for region III of Figure |^. These 
results support the simple interpretation that the concentrations of the dynamical species, 
the activator I~ and the inhibitor C102~, are key factors in the occurrence of a transverse 
instability. However, this picture may be overly simplified: for values of gel width larger 
than the upper critical value, the concentrations of I~ and C102~ remain well within the 
instability interval of the variable [MA]^;^ analysis while the stationary state remains stable. 



3. Discussion 



In this section, we have explored the parameter dependence of the Turing instability as a 
function of malonic acid boundary condition and gel width. The use of the CIMA reaction- 
diffusion system (as opposed to the CDIMA system) and its corresponding boundary species 
(iodide, rather than iodine, and chlorite, rather than chlorine) by experimenters makes di- 
rect comparison of our numerics with their results difficult. Nevertheless, our numerical 
simulations display similar features to those present in experimental results as control pa- 
rameters are varied. These are: (1) the transition to patterns with progressively larger 
numbers of fronts, and (2) the appearance and subsequent vanishing of a transverse insta- 
bility. Also agreeing with experimental results, the instability is localized near a minimum 
of the starch triiodide (as well as iodide) unperturbed solution. For the parameters and 
boundary conditions we considered in our numerical investigations, the multiply-peaked lin- 
ear instability eigenvectors remain stable, and the experimentally observed transition from 
a single to multiple unstable layers is not obtained. Overall, the agreement between trends 
in these experiments and our numerics is encouraging, and will hopefully provide motivation 
for future experiments on the CDIMA system which could be compared quantitatively with 
numerical and analytical results. 
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VI. CONCLUSIONS 



In this work, we have focused on the one-dimensional patterns that form in the presence 
of feed gradients, a necessary feature of the real experimental systems. We have shown 
that longitudinal diffusion along the boundary feed gradients can be important over length 
scales longer than the Turing wavelength. Therefore, the frequently-invoked locally uniform 
approach for predicting the linear instability of the stationary patterns along the gradients 
to a transverse symmetry-breaking instability is inappropriate in such cases. 

We have also explored the dependence of the Turing instability of these longitudinal 
structures on two control parameters. The transition to patterns with progressively larger 
number of longitudinal fronts and the appearance and subsequent vanishing of the transverse 
instability near a local minimum of the starch triiodide solution are features which are 
in agreement with experimental results. We have attempted to interpret these trends, by 
determining that a transverse instability occurs and is localized at that part of the stationary 
solution along the gradients where the values of the concentrations of the dynamical iodide 
and chlorite species are within a certain well-defined range. For the parameters and boundary 
conditions investigated here, we do not obtain the experimentally observed transition from 
a single to multiple Turing unstable layers. 

Building on the work presented here, one can begin to address many interesting ques- 
tions. At the linear level, it is clear from the numerical solutions that both the formation of 
localized structures along the gradients as well as their transverse instability are governed 
by the Turing mechanism. The relationship between the "longitudinal" and "transverse" 
pattern formation is an interesting question that so far has only been analyzed for model 
systems |^|. For the realistic chemical description it may be more feasible to address this 
question analytically in the limit where the reservoir concentrations can be approximated by 
linear profiles: we have found that these conditions are achieved, for parameter values ex- 
plored in the variable gel-width numerical investigation, for small gel widths. A description 
of the longitudinal structures can perhaps be sought, in terms of wave number selection in 



the presence of control parameter ramps The effect of the variation in background con- 
centration profiles, with changing gel width or malonic acid reservoir concentration (Figures 
p!0| and |l^), on these one-dimensional structures can be investigated within this framework. 
Alternatively the symmetry breaking transition might be better understood in terms of a 
fully three dimensional periodic instability perturbed by longitudinal gradients. 

At the nonlinear level, the nature of the transition from one-dimensional non-symmetry 
breaking front patterns to symmetry-breaking transverse spots in thin-strip reactors (con- 
tinuous or discontinuous) can be determined numerically from the nonlinear evolution of 
the LRE model in two dimensions Numerical establishment of the bifurcation behavior 
for the two-dimensional "monolayers" in disc reactors is also a relevant topic for further 
investigation, and would require extending numerical computation to three dimensions. Du- 
fiet et al. have pointed out that these monolayers which are confined by a transverse 
parameter ramp must be distinguished from what they call genuine two-dimensional struc- 
tures which form under uniform control parameters. They have compared pattern selection 
in genuine two-dimensional systems and in such monolayers in the context of an abstract 
react ion- diffusion model. The advantage of using the LRE model is that results would then 
be directly comparable with experiments based on the CDIMA reaction. An interesting 
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feature of the work of Dufiet et al. is the couphng of the pattern forming modes to the 
longitudinal displacement of the pattern as a whole: it would be interesting to investigate 
this effect using the realistic description of the longitudinal gradients. 

The successful formulation of a realistic model of the chlorine dioxide-iodine-malonic 
acid reaction-diffusion system has made this system an attractive paradigm for the study of 
nonequilibrium pattern formation [|l^]. This work represents an attempt at bringing theo- 
retical and numerical study of pattern formation in chemical systems closer to experimental 
studies. 
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APPENDIX A: DETAILS OF THE TURING INSTABILITY CONDITIONS 

1. Stability to Homogeneous Perturbations: k = 

The characteristic equation for A; = is: 

aXo^ - (an + cra22)Ao + (a22aii - 012021) = 0, (Al) 

with, 

Xo^^^ = 77-(aii + 0-^22) ± 77-\/ (cm + 0-^22)^ - 4o-(aiia22 - 012021)- (A2) 
We require Re(Ao*^^'') < 0. Therefore, in the case of Im(Ao*^^-') 7^ 0, we must have: 

Ao^+^ + Ao^"^ = an + aa22 < 0, (A3) 
and additionally, in the case of Im(Ao(±)) = 0: 

Ao^^'^Ao^^) = a [ana22 - 012021] > 0. (A4) 



2. Instability to Inhomogeneous Perturbations: A; / 

We require that at least one of the roots be positive for some A; 7^ 0. Consider the sum 
of the eigenvalues: 

Afc^+) + Afc(-) = -{ac+ l)k^ + (an + ^022). (A5) 

Once stability to A; = perturbations is imposed according to Eq. ( |A3D , the above sum will 
be < 0, and the real part of one of the roots is necessarily negative. First, we require that 
the product of the roots be < in order to have a positive real root: 
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Afc^^''Afc^~^ = a {ck^ - (022 + caii)P + (011022 - 0120-21)) = h{k'^) < 0. (A6) 
Since the first and third terms in h{k'^) are > 0, a necessary condition for h{k'^) < is: 

022 + can > 0. (A7) 
Furthermore, we know that has a minimum, and from Eq. ( [A^ ) that: 

h{0) = cr(aiia22 - 012021) > 0. (A8) 
So, for h^k"^) < 0, k'^ must he between the two roots k_^ and k^'^: 
1 



^ 2c 



(A9) 



(022 + can) ± (022 + caii)2 - 4c(aiia22 - 012021) 

For k±^ to be real, we must have: 

(022 + coii)^ - 4c(aiia22 - 0.12021) > 0. (AlO) 



Eq. ( |A10| ) is required for real values k± of the marginal modes, and Eq. ( |A7| ) is required 



for k±^ > 0. Equations (^), (^), (^) and ( |Aiq ) constitute the Turing conditions. 

3. Oscillatory Instability 

From Eq. (|^) there will be a complex conjugate pair of eigenvalues, Im(Afc) 7^ 0, for: 
^(P) = (ac - ifk^ + 2k^{ac - l)(an - (7022) + (on - (7022)^ + 4(Tai202i < 0. (All) 
To determine the behavior of gik"^), look at its roots: 



^(H)2 ^ (on - ff022) ^ V-4(TOi2a2i .^^2) 
^ (TC — 1 (7C — 1 ' 

For Im(Afc) 7^ 0, require fc^"* > 0. Under typical experimental conditions, we expect 
crc — 1 > 0. We note that glk"^) possesses a minimum with: 

g{0) = (an - (7022)^ + 4a-ai202i. (A13) 

Therefore, for fc^"* > 0, we must have g{0) < 0, in which case ImA^ 7^ for < A;^ < /c^^^ . 
The real part of the complex conjugate pair is given by: 

Re Afc = [{ac + l)k^ - (an + ^022)] , (A14) 

and behaves as: 

i) an + aa22 < =^ Re < (A15) 

tt) an + (ra22 > Re A^ < for fcf > k^ > ^ "^"'^ (A16) 

crc + 1 

ReA,>0 for < A:^ < ^il±^. (A17) 

- ac + 1 ^ ^ 
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RESERVOIR II 




FIG. 1. Sketch of open reactor geometries, adapted from Refs. [10, IJ]: A block of gel, with 
dimensions L ^ w > h is in contact with two reservoirs I and II. The reservoirs are continuously 
stirred and fed with fresh supplies of reactants, such that each is separately nonreactive. A gra- 
dient in the reservoir species forms perpendicular to the feed boundaries in the z-direction. The 
symmetry-breaking patterns form transverse to this gradient. The thickness A of the finite region 
along z where they form is equal to at least one wavelength A of the Turing patterns. 
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FIG. 2. One-dimensional steady state solutions: The boundary conditions are 

[MA]^ = 1 X 10~2 M at the left boundary, and [12]^ = 1 x lO"^ M and [C102]^ = 6 x 10"^ 
M at the right boundary. All other boundary conditions are zero. The z-axis has been normalized 
with respect to the thickness of the gel in the z-direction, w = 0.3 cm. We note that the steady 
state profiles for the reservoir variables malonic acid (MA), iodine (I2) and chlorine dioxide (CIO2) 
vary considerably from diffusion-only linear profiles. The series of peaks in the starch triiodide 
(8X3^) profile correspond to experimentally observed stripes parallel to the feed boundaries ||lF 
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FIG. 3. The value of the steady state solution, locally uniform in the x — y plane, is plotted at 
each point z for the two dynamical variables iodide and chlorite (C102~) using the reaction 
+ diffusion profiles of the reservoir variables given in Figure |2|. 
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FIG. 4. Boundaries from locally uniform stability analysis: At each location z, the Turing 
instability conditions have been plotted for the corresponding uniform steady state. The plotted 
quantities have different dimensions, and since only the sign of each quantity is of interest, the 
vertical axis has no scale. When the light solid line corresponding to (an + aa^-i) is less than 
zero, the complex conjugate pair of eigenvalues for A; = has a negative real part. The dotted 
line is {a\\a22 — 012021) which is additionally required to be greater than zero for stability of a 
real A; = mode. We note that this quantity is everywhere greater than zero. The long-dashed 
line corresponds to (can + 022), and the heavy solid line is (can + ^22)^ — 4c(aiia22 — 012021), 
both of which must be greater than zero in order to have a A; 7^ instability. The dashed line is 
((oii — cra22)^ + 4crai2a2i), and where it is less than zero, a complex conjugate pair of eigenvalues 
exists for a finite range in k. The shaded region indicates where the locally uniform steady state is 
stable to homogeneous perturbations and unstable to inhomogeneous perturbations. At z = 0.354, 
the light solid (an + C7a22) and long-dashed (can + 022) lines go through zero, while the heavy 
solid, dotted and dashed lines remain positive. 
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FIG. 5. Gain curves for focally uniform steady states: The solid line denotes the eigenvalue 
with the larger real part, A+, and the dashed line denotes the one with the smaller real part, A_. 
At (a) z = 0.20, both eigenvalues are real and negative for all k; (b) z = 0.50, both eigenvalues 
are real and peaked at A; = 0, with A_|- > in the range < /c^ < A;+^ and A_ > in the range 
< A;^ < kJ^ {k±^ are given in Appendix (c) z = 0.62, both eigenvalues are still real, but at 
A; = 0, we have a real degenerate pair, giving the boundary of the Hopf bifurcation; (d) z = 0.70, 
complex conjugate pair of eigenvalues with positive real part for < /c^ < (on + aa22) / {o'D^ + Du), 
and k ^ instability for < k"^ < k^^; (e) z = 0.72, same as (d), except that the real part of 
the complex conjugate pair is peaked at zero growth rate; (f) z = 0.74, same as (d) with the real 
part of the complex conjugate pair less than zero at A; = 0; (g) z = 0.76, same as (f), except that 
the maximum growth rate for A; 7^ is zero; (h) z = 0.80, same as (g), except that the maximum 
growth rate for A; 7^ is negative; (i) z = 0.85, both eigenvalues are real and negative. 
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FIG. 6. Gain curve for the one-dimensional steady state along the gradients: The real parts of 
the first fourteen eigenvalues with largest real parts have been plotted. The eigenvalues are real, 
except at points (or along intervals) where two curves intersect. The eigenvalue crossings appear 
imperfect due to the coarse selection of k values. All eigenvalues are less than zero, with the heavy 
line corresponding to the slowest decaying mode at each k. 



28 



(D 

■ I—' 

O 



CD 

■g 
o 




-8 X 10"5 



FIG. 7. Eigenvector corresponding to slowest decaying mode in a fully nonuniform analysis: 
The eigenvector at k = 81.6 cm~^ with the largest real eigenvalue has been plotted. 
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FIG. 8. Modified local stability analysis: K' , Hi and H2 are given by the solid, long-dashed 
and dotted lines respectively. In (a), these boundaries have been modified to take into account 
difi'usion of the steady state along whereas in (b) they are unchanged. The Turing instability 
region is indicated by the shading. The shaded region in (b) is identical to that in Figure ^ Note 
that the right boundary of the Turing region, which denotes ^ criticality, remains unchanged 
under the modification. 
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FIG. 9. Malonic acid concentration (M) at the left boundary versus the wave number k of the 
transverse instability (1/cm): The crosses represent the marginally stable wave numbers determined 
from the stability analysis of the one-dimensional steady states corresponding to the given malonic 
acid feed concentration. A transverse instability occurs for values of malonic acid feed concentration 
in the shaded regions, with the range of linearly unstable modes delimited by the solid lines for 
each value of [MA]^. 
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FIG. 10. Stationary solution for [Sl3~] along feed gradients for various values of malonic acid 
feed concentration at the left boundary: This series of plots shows how the stationary state along the 
gradients changes as [MAJ^ is varied. The number of peaks in the solution increases with increasing 
malonic acid concentration, while the left-most peak becomes smaller. The minima correspond to 
light bands in the experimental results. Instability region I corresponds to stationary states with 
three peaks, region II with four peaks and region III with five peaks. 
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FIG. 11. Density plots of the stationary solution along the gradients and the fastest-growing 
instability eigenvector for region II: White and black correspond to large and small values of the 
solution, respectively. The top plot shows the variation of the stationary state in region II with 
increasing [MA]^,. The bottom plot shows how the most linearly-unstable eigenvector is singly 
peaked and tracks the next-to-last minimum in the solution. (Note: It is difficult to discern the 
last minimum in the top plot.) The "staircase" structure is due to the discrete sampling of [MAj^. 



33 



4 X 10" 



^ 3x10" 




■g 
x 
o 

b 

03 
C 



O 




c 

T3 
O 



(c) 



2 X 10"2 
1 X 10"3 



3 X 1 0"3 
S 2 X 10"2 

1 X 1 0"2 - 



(b) 



n.^^ ° ° ° ° 

axEoaoxoao o o o o 

AAA/VW^AAAAA A A A A 




0.025 



0.050 



03 

■g 

o 



o 

CO 



o 
O 




■g 

o 



6 X 10"=^ 
4 X 10"^ 
2 X 10"^ 


4 X 10"^ 

2 X 10"^ 



3 X 10"'' 
2 X lO"'' 
1 X 10"'' 





(f) 



(e) 



A^ 



A^ 




o o o » 



^mxcoxo o o o o 



4 fb 



'-'i-manixj o o o o 

^^^A^WVl A A A A 







0.025 



Malonic Acid at Left Boundary (M) 



0.050 



FIG. 12. Values of stationary profiles of the CDIMA chemical species at successive minima 
of the stationary [8X3^] solution as a function of [MA]l: The open triangles, circles, squares and 
diamonds represent the second, third, fourth and fifth minimum, respectively. The filled symbols 
correspond to the critical values of [MA]l; the points between the filled symbols correspond to the 
linearly unstable states. We note that the instability occurs successively along the second, third 
and fourth minimum. Although the values of [MA], [I2], [CIO2] in the unstable ranges continue 
to increase as functions of [MA]/, in going from one instability region to the next, the trend for 
the dynamical species [I~] and [C102~] is different: the value of [C102~] remains approximately 
constant, while [I~] varies within an approximately constant range. 
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FIG. 13. Density plot of the iodide eigenvector corresponding to the largest eigenvalue as a 
function of the transverse wave number k and spatial distance along the gradients ([MAJ^ = 0.023 
M): Black and white correspond to low and high values of the eigenfunction, respectively. For small 
values of k, the eigenvector is multiply peaked while for k larger than approximately 300 cm~^, 
which includes the unstable range of wave numbers, it is singly peaked. 



35 



0.01 




k (cm) 

FIG. 14. Spectrum of eigenvalues for [MA]^ = 0.023 M: The real parts of the first fourteen 
eigenvalues with largest real parts have been plotted. The eigenvalues are real, except along 
intervals where two curves overlap. Some eigenvalue crossings appear imperfect due to the coarse 
selection of k values. It appears that the first two eigenspectra do not cross but remain distinct. 
We note that the eigenvalue crossings define distinct "modes" that cut across the spectrum of 
eigenvalues. The filled triangles indicate the eigenvalues for one such "mode" . 
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FIG. 15. Stationary solution for [Sis"] along feed gradients for various gel widths: This series 
of plots shows how the stationary state along the gradients changes as the gel width is varied. The 
horizontal axis is the scaled length along the gradients. With decreasing gel width, we note: (i) 
shifting of the pattern to the right, and (ii) increase in the concentration scale by approximately 
an order of magnitude (primarily due to the leftmost peak). 
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FIG. 16. Gain curve corresponding to largest eigenvalue for various gel widths: (b)-(d) are 
unstable; in (g)-(i), the largest eigenvalue now occurs at A; = 0. 
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FIG. 17. Gel width w (mm) as a function of wave number k of the transverse instabihty (1/cm): 
The crosses represent the marginally stable wave numbers determined from linear stability analysis 
of the one-dimensional steady states corresponding to the given malonic acid feed concentration. 
A transverse instability occurs for values of gel width in the shaded regions, with the range of 
linearly unstable modes delimited by the solid lines for each value of w. The vertical plot range 



corresponds to the experimental range of gel width in the bevelled thin-strip reactor |14]. 
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FIG. 18. Values of stationary profiles of the CDIMA chemical species at successive minima of 
the stationary [Sl3~] solution as a function of gel width, w: The open triangles, circles, squares and 
diamonds represent the second, third, fourth and fifth minimum, respectively. The filled symbols 
correspond to the critical values of gel width; the points between the filled symbols correspond to 
the linearly unstable states. 
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TABLE I. Kinetic constants for the CDIMA system. 



Rate or diffusion constant 


Dimensions 


Value 


ha 


(s-i) 


9 X 10"'' 1 


hb 


(M) 


5 X 10^^ ^ 


h 


(M-^s-i) 


1 X 10^ ^ 


ha 


(M-2s"i) 


1.2 X 10^ ^ 


hb 


(s-i) 


1.5 X 10"^ ^ 


h 


(M2) 


1.0 X 10"^^ ^ 




(M-2g-l) 


6.0 X 10^ 2 




(s-i) 


1.0 2 


Di- 


(cm^s"^) 


7.0 X 10~6 3 


-^ClOa" 


/ 9 _1 \ 

(cm s ) 


7.0 X 10~^ ^ 




(cm^s""*") 


6.0 X 10^^ 1 


Dma 


(cm^s""*") 


4.0 X 10"^ 1 


Dc\02 


(cm^s^"*") 


7.5 X 10"^ 1 


Dh+ 


(cm^s"^) 


1.0 X 10-5 


K[S]o 


(M-i) 


6.25 X lO'' ^ 


iProm m at 7°C; ^From |7 


at 4°C; ^From [|| at 4°C; ^From [| 


1 at 4°C. 
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